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(-^ \ A nonlinear theory of pattern selection in parametric surface waves (Faraday waves) is 

presented that is not restricted to small viscous dissipation. By using a multiple scale 
asymptotic expansion near threshold, a standing wave amplitude equation is derived 
from the governing equations. The amplitude equation is of gradient form, and the 
coefficients of the associated Lyapunov function are computed for regular patterns of 
various symmetries as a function of a viscous damping parameter 7. For 7^1, the 
selected wave pattern comprises a single standing wave (stripe pattern). For 7 C 1, 
patterns of square symmetry are obtained in the capillary regime (large frequencies). At 
lower frequencies (the mixed gravity-capillary regime) , a sequence of six- fold (hexagonal) , 
eight-fold, . . . patterns are predicted. For even lower frequencies (gravity waves) a stripe 
pattern is again selected. Our predictions of the stability regions of the various patterns 
are in quantitative agreement with recent experiments conducted in large aspect ratio 
c/2 I systems. 
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1. Introduction 



This paper extends an earlier calculation by [Zhang fc Vinals (1997)| of the amplitude 
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equation governing Faraday waves in the weakly nonlinear regime. In order to make 
the problem analytically tractable, they neglected without rigorous justification viscous 
terms in the boundary conditions at the free fluid surface that had a nonlinear dependence 
on either the surface displacement away from planarity, or on the surface velocity. Even 
though the resulting amplitude equation led to the prediction of stationary patterns that 
are generally in agreement with experiments conducted in the regime of weak viscous 



dissipation ( Kudrolli fc Gollub (1996) Binks fc van de Water (1997) ), the unsystematic 



nature of the truncation makes it difficult to asses the range of validity of the theory. In 
particular, the so-called stripe pattern (a pattern comprised of a single standing wave) 
which is generically observed when viscous dissipation is not small, could not be obtained 
in their analysis for any range of parameters. We extend below this earlier work, and 
present a systematic weakly nonlinear theory of Faraday waves. Our results on pattern 



selection agree with those of Zhang & Vinals (1996) and (1997) in the limit of small 



viscous dissipation, and with recent experimental work otherwise. 

Parametrically driven surface waves (also known as Faraday waves) can be excited on 
the free surface of a fluid layer that is periodically vibrated in the direction normal to the 
surface at rest when the amplitude of the driving acceleration is large enough to overcome 



the dissipative effect of fluid viscosity (Faraday (1831), Miles & Henderson (1990)). Of 
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special concern to us is the issue of pattern selection in a layer of lateral dimension much 
larger than the excited wavelength (see, e.g., Cross & Hohcnbcrg (1993) for a recent re- 
view on pattern formation). In the case of Faraday waves, it is now known that different 
wave patterns can be excited depending on the fluid properties and the driving ampli- 
tude or frequency. At high viscous dissipation (a fluid of large viscosity and/or a low 
driving frequency) , the observed wave pattern above threshold consists of parallel stripes 
(Edwards & Fauve (1994) and Daudet et al. (1995)| ). For lower dissipation, patterns of 

square symmetry (combinations of two perpendicular plane waves) are observed in the 

capillary regime (large frequencies) (Lang (1962),[Ezcrskii, Rabinovich, Rcutov fc Starobincts (1986) 






Tufillaro, Ramshankar fc Gollub (1989)1 pilibcrto, Douday fc Fauve (1991) 



Christiansen, Alstr0m fc Lcvinscn [1992] , |Miiller (1993)| , [Edwards fc Fauve (1994)[) . At 



lower frequencies (the mixed gravity-capill ary regime), hig her symmetry patterns have 
been observed by Kudrolli fc Gollub (1996) (hexagonal) and Binks fc van dc Water (1997) 
(hexagonal, eight- and ten- fold) . The aim of this paper is to present a weakly nonlinear 
analysis of Faraday waves that predicts stationary wave patterns with these symmetries. 
The derivation of an amplitude equation is a classical method to describe excited states 
beyond linear instability. Just above threshold, the evolution of the system is assumed to 
be described in terms of the complex amplitude A of the most unstable mode according 
to linear theory. The equation of motion for A is often of the form 

^=*A-gA\ (1.1) 

where a is the linear growth rate, and g > is real. The low order nonlinear term 
provides saturation. There exist cases, however, in which spatial isotropy permits waves 
to be excited in any direction, and the nonlinear interaction term in the equation above 
contains terms of the form gij\Aj\ 2 Ai, with Ai and Aj the slowly varying amplitudes of 
two degenerate unstable modes. If the coupling coefficients gij are known, the resulting 



wave pattern can be predicted from Eq. (1.1), as has been illustrated by Muller (1994) 
for Faraday waves. 

The derivation of amplitude equations for surface waves is greatly simplified in the case 
of an ideal (inviscid) fluid. Since the bulk flow is irrotational, there exists a hamiltonian 
formulation in which the canonically conjugate variables are the surface displacement 
and the velocity potential at the free surface ( Zakharov (1968)] Miles (1977)). As a 
consequence, early analyses of Faraday waves were based on the hamiltonian descrip- 
tion of the inviscid limit, and treated viscous or dissipative effects as a perturbation 
flMiles (1984)1 , |Milner (1991)], |Milcs (1993)D- The derivation usually starts from the set 
of ideal fluid equations (Benjamin & Ursell (1954)), written in terms of the surface veloc- 
ity potential 4>- The linear or zero-th order solution 0o is a sum over waves of frequency 
u> and wavevector {kj}, with u) and k = ||kj || related by the ideal fluid dispersion relation 
(Eq. (pT3T ) below): 



~kz 



k 



J2^(T)e i{kr ^ UJt) +- 



where T — et is a slow time scale, with t C 1 the dimensionless distance away from 
threshold. An expansion o f the ideal fluid equations to third order in e yields the equation 
for Aj(T) QMilner (1991)| ) 



dAj ikf 






dT 



(1.2) 



with / the amplitude of the driving acceleration, and II real functions of the angle between 
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the j'-th and k-th wavevectors. Since the coefficients of the cubic terms are imaginary, 
these terms do not contribute to wave saturation. This can be seen, for example, by 
computing d\Aj\ 2 /dT and observing that all cubic terms cancel. This is also a man- 
ifestation of a general symmetry principle in hamiltonian (or reversible) systems that 
prohibits real coefficients of cubic nonlinear terms in standing wave amplitude equations 
(pross fe Hohenberg (1993], poullet, Frisch fc Sonino (1994|). 



In the limit of low viscous dissipation, Hamilton's equations have been supplemented 
with a dissipation function ( [Miles (198TJ| , [Milner (1991)1 [Miles (1993)] ), which is com- 
puted under the assumption that the dominant contribution to viscous dissipation arises 
from the irrotational velocity field in the bulk, and not from friction at the container walls 
or dissipation near the free surface (where vorticity is produced) . Under this assumption, 
the rate of energy loss is given by (Landau & Lifshitz (1959)), 



E = -2r, dV 



U 2 



dxidxj 



where 77 is the shear viscosity, and the integral extends over the bulk fluid. The velocity 
potential <f> is now expanded in powers of e, and viscous contributions computed order 
by order in e. This procedure leads to imaginary components in the coefficients of the 



cubic terms of Eq. (1.2), and therefore to wave saturation. The precise functional form 



of the coefficients o btained by this method is still somewhat controversial (Miles (1993) 
Hansen fe Alstr0rr| ). 



We next address the effect of the rotational component of the flow. The dimension- 
less group involving the ratio of viscous to inertial effects is the damping parameter 
7 = 2^fc 2 /cjo, where fco is the critical wavenumber in the inviscid limit, and ujq its 
angular frequency (7 is inversely proportional to the Reynolds number of the flow). 



Lundgren & Mansour (1988) considered an expansion of the governing equations and 



boundary conditions in powers of 7. They showed that in the weak dissipation limit, the 
dominant terms in the boundary conditions are 0(7), with a first correction at 0(j 3 ^ 2 ). 
At linear order in the surface displacement or surface velocity, terms of 0(7) are purely 
irrotational, while the rotational flow component contributes at 0(7 3 / 2 ). In fact, linear 
stability analysis of Faraday waves by Miiller et al. (1997) (see also Sec. pj) shows that 
the dimensionless value of the driving amplitude at threshold equals 7, with a first correc- 
tion term that is proportional to 7 s ' 2 . The dominant contribution arises solely from the 
irrotational flow com ponent, with contributio ns from the rotational component coming 
at C(7 3 / 2 ). However, Zhang fc Vinals (1997) argued that the lowest order contribution 
from both irrotational and rotational components is 0(7) at the nonlinear level in the 
surface variables. Hence rotational flow cannot be neglected in a nonlinear theory, even 
in the limit of small dissipation. For example, the kinematic boundary condition at the 
free surface does include one such term proportional to 7 that arises from the component 
normal to the surface of the rotational part of th e velocity field. This term was retained 
both in the analysis of Zhang fc Vinals (1997)| and in our analysis below, but not in 
previous approaches based on a dissipation function. 

Another qualitative feature of importance to pattern selection in Faraday waves is 
triad resonant interactions. Since the standing wave amplitude equation must be invari- 
ant under Aj — > — Aj, m triad resonance cannot contribute directly to quadratic order in 



| The governing equations are invariant under time translation by a period of the driving 
force t — » T + 2n/ui. Subharmonic response implies that ((x,y,t + 2tt/lu) — —£(x,y,t), with 
z = £(x, y, i) the position of the surface. Because of this invariance, the amplitude equation for 
Aj must also be invariant under a sign change in Aj 
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Aj , but it does contribute significantly to the coefficients of the cubic order term via the 
coupling of the zero-th order unstable wave and first order stable waves. This resonance 
was already encountered by Milner (1991)| as a divergence of the cubic coefficient in the 



standing wave amplitude equation at a particular angle. Later, Edwards & Fauve (1994) 



suggested that triad resonance would be important at low viscous dissipation, range in 



which linearly stable modes are only weakly damped. Zhang & Vinals (1997) calculated 
such a contribution explicitly, and showed that it is important in determining the sym- 
metry of the selected pattern in the region of small 7. In particular, they predicted a 
sequence of quasi-periodic patterns in the region in which the resonant angle approaches 
zero. 



We extend in this paper the analysis of Zhang & Vinals (1997) that was based on a 



quasi-potential approximation to the governing equations. By separating the rotational 
flow within a small vortical layer near the free surface from the potential flow in the bulk, 
they derived a standing wave amplitude equation valid in the limit of small viscous dis- 
sipation. The calculation, however, relied on an uncontrolled approximation concerning 
nonlinear viscous terms and, as a consequence, its region of validity is difficult to asses. 
We describe below a systematic expansion of the Navier-Stokes equation and boundary 
conditions that overcomes this difficulty and that leads to an amplitude equation not 
restricted to small viscous dissipation. We start by deriving an exact, although implicit, 
relation for the threshold of instability, which is then used in the nonlinear analysis. This 



result extends earlier numerical work by Kumar fc Tuckerman (1994) , and agrees with 



a recent low viscosity approximation to the governing equations by Midler et al. (1997) 
We then use a multiple scale expansion to derive a standing wave equation which is of 
gradient form. Minimization of the associated Lyapunov function leads to the predic- 
tion of stationary patterns of different symmetries as a function of the fluid parameters 
and frequency of the driving acceleration. Our predictions are in good agreement with 
experiments conducted in large aspect ratio cells. 



2. Governing equations and linear stability 

We consider a semi-infinite fluid layer, unbounded in the x — y direction, extending 
to z = —00, and with a planar free surface at z = when at rest. The fluid is assumed 
incompressible and Newtonian. Under periodic vibration of the layer in the direction 
normal to the surface at rest, the equation governing fluid motion (in the co- moving 
reference frame) is 

d t u+(u- V)u= — Vp + is\7 2 u + g z {t)e z , (2.1) 

P 

with u the velocity field, p the pressure, p and v the density and kinematic viscosity of 
the fluid respectively, and g z (t) = —g — fcosuit the effective gravity, [j] The base state 
is a quiescent fluid with a pressure distribution p = pg z (t)z. We absorb the body force 
in the pressure, so that in what follows p is the deviation from pg z {t)z. By applying 



— (V x Vx) to Eq. (2.1), one can eliminate the pressure term, and also obtain a system 
of equations for the velocity components of u = (u,v,w), in which the linear terms are 
uncoupled, 

<9 t V 2 u - ;A7 2 V 2 u = V x V x (u ■ V)u. (2.2) 

t We use a drivi ng acceleration proportional to cos uot instead of sin uit as in 



Zhang & Vinals (1997) to avoid, as discussed in that reference, the complication related to 
the parity under time reversal of the driving acceleration. 



Amplitude equations and pattern selection .. 



Continuity, V • u = 0, has also been used to derive Eq. (2.2). 

Besides the null conditions at z — — oo, there are four boundary conditions at the 
moving free surface (Lamb (1945)). Let z = ((x,y) be the position of the surface (Fig. 
0), then the outward pointing unit normal fi is 



(-g,C, -dyS, 1) 

[i + (cU) 2 + (W] 1/2 ' 

Two linearly independent tangential unit vectors ti and £2 are 

t _ (1,0,9,0 t _ (°>W) 



[l + (9xC) 2 ] 1/2 ' 



[i + ao 2 ] 1/2 ' 



(2.3) 



(2.4) 



Note that these two vectors are not mutually orthogonal. The choice is made so that the 
expressions for ti and £2 are symmetric in the Cartesian variables x and y. 
The kinematic boundary condition is, 

d t (+(u(z = 0-V H )( = w(z = (), 

with Vjj = e x d x + e y d y . Since the governing equations will be expanded and solved 
order by order, we quote here its Taylor expansion around z = 0, 



d t ( +[u + d z uQ z=0 d x ( +[v + d z v(;} z= ody( = [w + d z w( + ^d zz w( 2 } z=0 



(2.5) 



Only terms up to third order in the velocity or surface displacement will be required, f 
Neglecting the effect of the air phase above the fluid, the tangential stress at the free 
surface is zero, 

t m • T • n\ z=c = 0, to =1,2, 

with T the stress tensor of components, Tij — [—p — pg z (t)z]5ij + pv(djU{ + diUj). The 
normal stress at the fluid surface is balanced by capillarity, 

h-T-h\ z=c = 2Ha, 

where a is the surface tension and 2H is the mean curvature of the surface, 

2H = {d XX C [1 + (dyC) 2 ] + dyyC [l + (d^] - 2d X C,dy(,d X y(\ / [l + (d^ + (dyC) 2 ]^ ■ 

The li near stability of the fluid la yer under vibr ation was first addressed in the inviscid 



limit by Benjamin fc Ursell (1954', and later by Landau fc Lifshitz (1976) in the limit 
of low viscosity. More recently, Kumar fc Tuckerman (1994) numerically compute d the 
neutral stability curve for a fluid of arbitrary viscosity, and Muller et al. (1997) have 
given an analytic low viscos ity expansion. We first review briefly the formulation of 
Kumar fc Tuckerman (1994) , and then show that an exact (albeit implicit) analytical 
expression for the threshold can be derived, thus avoiding the numerical calculation. 

The domin ant response of the param etrically driven system is sub-harmonic at a fre- 
quency w/2 ( Benjamin fc Ursell (1954) ). Although the methodology discussed below 
can also be used to analyze a possible harmonic response, we restrict our analysis to the 
subharmonic case. To address the linear stability of the fluid surface, we consider the 
following solutions for the vertical velocity field and surface displacement, 



w = cos(£;:e) ^2 e jiuJt/2 w J (z)Aj + c.c. 

j'=l,3,5,~ 



(2.6) 



f In order to avoid excessive use of parentheses, we follow the convention in the remainder of 
the paper that the operator d acts only on the function immediately following it. 



Co = cos(kx) 2_] 

3 = 1,3,5 
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e liut/2 A +c.c. 



where the Aj are complex amplitudes, and we retain all the harmonics of the fundamental 
mode e 1 "'/ 2 . Truncation of the sums to include the fundamental mode e luJt l 2 alone is 



only appropriate for small viscous damping. From Eq. (2.2), the linearized equation of 
motion for u>q is 

(<9 t V 2 - ^V 2 V 2 ) w a = 0. 

Substituting wq and Co from ( |2.6|) , one finds 



\jiu{-k 2 + 8 ZZ ) - v (-k 2 + 8 ZZ ) 2 \ w J (z) = 0. 

The solution of this equation is a linear combination of e ±kz and e ±qjZ , with q 2 = 
k 2 + jiui/2v. The linearized kinematic and tangential stress boundary conditions are 



dtCo 
(V 2 , - d z 



w 



By using the boundary conditions (pj) and the null conditions at z = — oo, Wq(z) is 
given by 

tuj»(z) = \-jiu + 2vk 2 ) e kz - 2vk 2 e q * z . 



The first term on the right hand side is the irrotational component of the flow, in which 
we have explicitly separated the inviscid and viscous contributions. The second term in 
the right hand side is the ro tational compo nent (this is the component that has been 
neglected in earlier work by Milner (1991) and Miles (1993) ). The linearized normal 
stress boundary condition is, after having eliminated the pressure by using the equation 
of motion. 



[2*A7 2 



H 



(d t - i/V 2 )] d z w + (g 



-Vn + fcosuit 

p 



VKo = o. 



(2.8) 



By substituting the assumed solutions given by Eqs. ( |2.6| ) into this equation, we note 
that the term 2vV 2 H d z WQ when acting on the irrotational flow component e hz yields a 
contribution at low viscosity that scales as v, whereas the rotational contribution (from 
e qjZ ) scales as v^l 2 . The remaining term {p t — ^V 2 ) d z Wo is simply equal to — uj 2 . Hence 
it is justified to neglect the rotational flow component in the linear stability analysis at low 



damping. As we show below, and in agreement with the work by Miiller et al. (1997) 



rotational flow contributes terms of order i/ 3 / 2 and higher to the value of the driving 
acceleration at onset. In the analysis that follows, however, we will retain the full linear 
solution. 

Kumar fc Tucker man (1994)| 



By equating the coefficients of each harmonic e^ luJt ^ 2 resulting from Eq. (2 
found 

H 1 A 1 - fA\ - fA 3 = 0, 
H 3 A 3 - fA 1 - fA 5 = 0, 
H 5 A 5 - fA 3 -fA 7 = 0, . . . 



(2.9) 



with 



Hi 



[4^-fc 4 - k(q 2 



k 



2\21 



gk 2 



This is a system of equations in the unknowns Aj , function of wavenumber k and driving 
amplitude /. By truncating the system (2.9) at some particular A n , it can be solved 
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numerically as an eigenvalue problem, / being the eigenvalue. This is indeed what was 
done by Kumar & Tuckerman (1994). However we observe that after truncation at A n , 

fA n -4 



A, 



fA n - 2 



A n -2 = 



■"« Hn-2 

so that the infinite set of equations can be re-written as 



J2_, 



H t 



r- 



Ho 



H B — 



A 1 -fAt = H 1 (k,f)A 1 -fAt=0. 



For a given wavenumber k, its threshold of instability /& is given implicitly by 



fh = \Hi(Jt,f k )\ 



#1 - 



fl 



H* 



(2.10) 



(2.11) 



(2.12) 



The complex amplitude Aj can then be recursively obtained from Eq. ( 2.10 ) and ( 2.11 ) 
up to a real factor. The critical wavenumber for instability fc onsct corresponds to the 
lowest value of /&, Jo- 



lt is interesting to consider the limiting behavior of Eq. ( 2.12 ) at low viscosity. First 
r ecall that for a semi-infin ite inviscid fluid, the dispersion relation for surface waves is 
QLandau fc Lifshitz (1959)1 ) 

ul = gk + o-kl/ p, (2.13) 

with wo = w/2, and kg the wavenumber. In a fluid of low viscosity we expect fc nsct to be 
near kg. It is then convenient to define dimensionless variables by using 1/wo as the time 
scale, and 1/fco as the length scale. We also define a reduced wave number k = k/ko, 
a viscous damping coefficient 7 = 2vk'^/ujQ ) the gravity wave G = gk^/oj^ and capillary 
wave £ = ak\j puj^ contributions to the dispersion relation, and the dimensionless am- 
plitude of the driving acceleration A = fko/4uiQ. Note that G + £ = 1 from Eq. ( 2.13 ); 
G = 1 corresponds to a pure gravity wave while G — to a pure capillary wave. For 
7 <C 1, fc onsct and A onsct in Eq. ( 2.12 ) can be expanded as a power series of the damping 
coefficient 7 

1 1 2_,yV2 4. -7+2G 2 , 

1+ 3-2G' + (3-2G) 2 ' T-.., 

5/2 , 



"'OllSCt 

^onsct 



7 



h m 



JJ-2G_ „ 

8(3-2G) ' 



(2.14) 



The first correction term is proportional to 7 3 ' 2 , and agrees with a low viscosity expansion 
of the linearized equations given by Muller et al. (1997) As an example, we plot in Fig. 
H the value of the threshold, A on set, as a function of 7 at G = 1/3. Previous low damp ing 
calculations of the stand ing wave amplitude equation by Milncr (1991) , Miles (1993) and 
Zhang & Vinals (1997) considered the dominant term A onsot = 7 only. Note, however, 
that the first correction, — 17 3 ^ 2 , can be a sizable contribution even for small 7 (e.g., a 
15% difference at 7 = 0.1). As a reference, we note that a similar linear analysis based 
on an inviscid formulation to which viscosity is added through a dissipation function, 
leads to the the damped Mathieu equation, 

d^Ck(t) +ld t Ck(t) + ^ 2 (l + 2Acos2wo*)C*(*) = 0, 

where Cfe(^) is the Fourier transform of £(x, i). This equation gives a threshold at 7 + 
37 2 /64 + 0(-f 3 ), which is plotted as the dot-dashed line in Fig. S. The first correction 
term is of a different order and has a different sign. Finally, we mention that rotational 
flow at the linear level in the surface variables can be incorporated into the damped 



Mathieu equation, as shown by Nam Hong (1993) 



8 P. Chen and J. Vinals 

3. Standing wave amplitude Equation 



In this section, we use the multiple scale approach of Newell & Whitehead (1969) 



to 

derive standing wave amplitude equations valid near threshold. It is interesting to note 
that the solvability condition in this case arises from the boundary conditions, unlike most 
other problems. For a driving amplitude / above threshold, we define e = (/ — /o)//o 
and expand the flow as 



1/2 

e ' u 



eui 



3/2 



and similarly for p and C- Near threshold, i.e., for e < 1, we separate fast and slow time 
scales: T — et; d t — + d t + ecV- Spatial slow scales are not included because only regular 
patterns are considered here. At order e 1 / 2 we recover the linear problem discussed 
in the previous section. Since we are interested in standing wave patterns of different 
symmetries, the solution at this order is written as a linear combination of waves with 
wavevectors k m of magnitude fc onsot but along different directions on the x-y plane, 



(/!() 



= Y>os(k m -x)B ro (T) 



Co = ^2 C0S ( k ™ ' x )^m(T) 



E 

3=1,3,5, 

E 

3=1,3,5,- 



jiut/2 j / \ 



c.c. 



e> iwt/2 ei 



c.c. 



where B m (T) are real wave amplitudes, functions only of the slow time scale T, and the 
ej are the same as the Aj found in Eqs. (2.10) and (2.11). 
At order e the equation of motion is 



(d t V 2 - i/V 2 V 2 ) wi = [V x V x (u • V)u ] • e z . 
By using the linear solution, the right hand side of Eq. t pJ\ ) is of the form, 

e^%(z). 



in n 



COS ((km ± k„) • x) ^ 
J=0,l,2, 



(3.1) 



(3.2) 



The first summation is over all possible k m ± k n , except zero, and the hj(z) are com- 
binations of exponential functions. Since every term in the right hand side contains a 
perio dic function of x, and exponential functions of t and z, the particular solution of 
Eq. (3T), Wip, can be easily found. The homogeneous solution wih and £i 



Wlh 



^2 cos((k m ±k„) -x) 



E 

ntn 



cos((k m ±k„) -x) 



E 


eJ iut ^|k m± k„| Wm ± + e r&*pi±J + c. c . 


Li=l,2,3,. 






, |k m ±k„|« 0± , Jk m ±k„| 2/ oO± 
'^ "-run ~ ^° t-'mn 


E 


e "mn ' c,c - ' ®mn 


(3.: 


L .7=1,2,3,- 







must now satisfy the boundary conditions. We have defined (r^ n ) = |k TO ±k„| 2 H- i 7«w/V. 
The constants a mn ,(3 mn and S mn are determined by the boundary conditions. At this 
order the boundary conditions are, 

dtCi -Wi = Gii(u ,Co) 
(V^ -d zz )w! = Gi 2 (u ,Co) 

(~d t + 2,vV 2 H + vd zz )d zWl + (g - ~V| + Jo cosw^VffCi = Gi 3 (u , Co)- 

P 
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where the functions Gn, G12 and G13 are listed in appendix \A\ For each wavevector and 
harmonic (each m, n, and j in Eq. (3.3)), the three boundary conditions are sufficient to 
determine the three unknowns a mn , Pmn and S mn in ( |3.3[ ) . Because the algebra is quite 
involved (the number of terms is on the order of several thousand) , we have in practice 
used a symbolic manipulation package to solve for these constants. 
At order e 3 / 2 the equation of motion becomes, 

(d t V 2 - vV 2 V 2 ) w 2 = ~d T V 2 w Q + {V x V x [(u • V)m + (ui • V)u ]} • e 2 (3.4) 

At this order we only need to consider resonant terms; for example, all terms proportional 
to cos(ki • x). The right hand side of Eq.(|3.4|) is of the form, 



cos(ki-x) J2 e jiujt/2 Ej(z). (3.5) 

j=l,3,5,— 

where we have used the solutions (uo, Co, u i, Ci) already determined. Again, the Ej(z) 
are combinations of exponential functions. The solutions for u>2 and £2 are 



u>2 = cos(ki • x) J2 e jiut ^[E j (z) + (a j e kz + b j e q " z )C j 

3=1,3,5,- 

C2 = cos(k! • x) J2 e 3ll " t/2 Cj 

.3=1,3,5," 

Here Ej (z) is the particular solution that corresponds to the right hand side at this order 



shown in Eq.(3.5), and a 3 e fcz + bje qjZ is the homogeneous solution, which has the same 
form as the linear solution. We now use the kinematic and tangential stress boundary 
conditions at this order to determine the constants aj and bj , so that the normal stress 
boundary condition yields a solvability condition for the amplitudes Cj, which in turn 
leads to the amplitude equations for the B m . (Note that there are various terms of B m 

The boundary conditions at this order are 

d t ( 2 -w 2 =G 2 i(u ,Co,Ui,Ci) 
(y 2 H ~ dzz) w 2 = G 2 2(u , Co, ui, Ci) 

(-d t + 3v\7 2 H + vd zz )d z w 2 + (g- -V 2 H + /o cos ut)V 2 H ( 2 - G 23 (u , Co, Ui, Ci), 

P 

where the functions G21, G22 and G23 are listed in appendix |A|. By using the first two 
equations, aj and bj are found to be (again with the help of a symbolic manipulation 
package) 

a 3 C 3 = v(k 2 + qj)Cj + EJ 
b j C j = -2vk 2 C j +E b j . 

Here E? and E b - are complicated expressions involving the amplitudes of the waves, B m . 
Now w 2 and C2 are substituted into the third boundary equation to yield 



H X C X - 


- f c; - 


- /0G3 = 


= ^1 


H3C3 - 


- /<A - 


- /0G5 = 


= ^ 3 


H5C5 - 


- /0G3 - 


- /0G7 = 


= ^5 



The left hand side of this system of equations is identical to Eq. ( p.9[ ) for the linear 
problem, and the functions Fj on the right hand side are functions of B m and dBi/dT. 
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Solving for Cj just like we solved for the linear threshold, we obtain 

#iCi - f Cl =F 1 + ^-(f 3 + ^-(F 5 + -- .)\ = F, 



with Hj defined similarly to Hi in Eq. (2.11). Since the threshold of linear instability 



given by /o = |J?i|, a nontrivial solution for C\ will exist if the following solvability 
condition is satisfied: 

FH*+F*f Q = 0. 

This condition immediately leads to a standing wave amplitude equation for B\, 

dFl 

-^ r =aB 1 - g Bf - £ g(0 ml )B^B lt (3.6) 

with 9 m i the angle between k m and ki. The linear coefficient a (times e) is the linear 
growth or decay rate of this wave, and can be obtained from the linear analysis (simply 
consider an extra factor e at in Eq. (|2.6|)). The coefficient g{9) describes the nonlinear 
interaction between different linearly unstable modes, and provides for the saturation 
of the wave amplitude. Figures || and [| show our results for different values of 7 and 
for S = (pure gravity waves), and £ = 1/3 (mixed gravity-capillary waves). It is also 
important to note the asymptotic behavior of g(9) as v — > 0. We have already discussed in 
Sec. y, that the irrotational component of the flow contributes to order v to the linearized 
equation of motion (Eq. J2.S|) ), whereas the rotational flow contribution scales as v 3 ! 2 
instead. This observation is the basis for earlier low viscosity approximations in which 



only viscous d i ssipation aris ing from the irrotational flow was considered ( Miles (1984) , 



Milncr (1991)[ [Miles (1993)|) . However, we have computed the coefficient g{9) with and 
without the linear rotational flow and observed that both contributions are of order v at 
small v. Therefore, a formulation that does not incorporate the rotational flow explicitly 
cannot obtain the correct form of the third order damping coefficients, even in the limit 
of small viscosity. 

Particular nonlinear interaction terms that contribute to g(9) are shown in Fig. M Two 
linearly unstable modes with wave vectors k m and k„ (|k m | = |k„| = fc on set) interact 
to produce a wave at k m + k„ with an amplitude proportional to B m B n . This mode 
corresponds to a first order solution (wi and Ci in Eqs. ( fT^ ) and Eq. fl3.3|)). Now 
k m + k„ couples back to the original wave at — k n to give a contribution B 2 B m to 
dBm/dJ '. Since the mode k m + k„ is damped (only waves with wavenumbers near fc on set 
are unstable), this is a dissipative term and contributes to nonlinear saturation of the 
wave. Triad resonance occurs when the frequency of the mode k m + k„ equals the driving 
frequency (the modes k m and k„ oscillate at half the driving frequency). Energy is now 
directly transferred into this mode which can have a very large amplitude at low damping. 
Since k m + k„ couples back to — k„, it provides a dissipation channel for the mode k„. 
Dissipation is enhanced by triad resonance and results in a large value of g(9 mn ) in the 
vicinity of the resonant angle. The resonant angle can be estimated from the inviscid 
dispersion relation (2.13D, written in dimensionless form, 



^ = Gfc + Sfc d , (3.7) 

with k = 1, and u> 2 — G + S = 1 for the linearly unstable mode. At resonance, we have 
ui = 2, and the resonant wave number k r = |k m + k n | satisfies k r (G + T,k 2 j = 4. If 9 r 
is the resonant angle between k m and k„, k r = y/2{\ + cos 6V), the resonance condition 
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becomes 

V / 2(l+cos6» r )[G + 2(l + cos(9 r )£] =4. (3.8) 

Because G + S = 1, this condition can only be satisfied when X > 1/3. For example, 
cos6» r = 2 1 / 3 - 1 for S = l. 

For finite damping, the resonance condition is modified. However, triad resonance is 
expected to be significant only at low damping because of the damped nature of the first 
order wave. For example, Fig. O shows g(9) for different 7 and E = 1. At small 7, the 
nonlinear coefficient grows near resonance and peaks at the resonant angle. The value 
of the peak is seen to decrease with increasing 7. At 7 = 0.1, resonance has almost 
disappeared. 

In the formulation presented earlier, resonance arises from the homogeneous solutions 
wih and £1, which require finding the constants a mn ,(3mn and S mn in Eq. ( [3.3] ) by 
enforcing the boundary conditions at first order. The boundary conditions give rise to a 
system of linear equations for a mn , (3 mn and 5 mn , the left hand side of which (its matrix 
form is explicitly given in appendix pi) at 7 = has a determinant 

8fc 2 (Gk 2 + Efc 4 ) [4k - (Gk 2 + Sfc 4 )] 2 , 
which, when equated to zero, is equivalent to Eq. (13. 



4. Pattern selection and comparison with experiments 

Since the standing wave amplitude equation (|3.6[) can be written in gradient form, 



the selecte d p attern near threshold immediately follows (Cross & Hohenberg (1993)) 
Equation ( |3.6|) is equivalent to 

dB n ST 

~~dT~ ~~~5B^' 
with the Lyapunov function T given by 



^=-J«E B ™ + lEE 9{Vmn)B 2 n B. 



m n' 



with go = g{0 n n) which equals half the value of g(9 — * 0). The amplitude equation then 
implies that 

dT _ y-^ ST dB n _ y^ / dB n \ 

dT ~ ^ 6B n dT ~ 2^\dT J ~ ' 

so that the preferred pattern can be determined by minimization of T. The experimen- 
tally observed regular patterns above onset consist of N standing waves, with uniform 
amplitudes and wavevectors k m ,m = 1...7V. The case N = 1 corresponds to a single 
standing waves (a pattern of parallel stripes), N — 2 to a pattern of square symme- 
try, N = 3 of hexagonal symmetry, etc. For these regular patterns, the standing wave 
amplitudes are 

B » = . ^ a ur~\> n = l---N. 

The value of the Lyapunov function as a function of N then becomes 



,2 



N 



4 .9o + £ m=2 g(6> m i) 
Figure shows the computed values of J-(N) as a function of 7 for different values of 
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E. For pure gravity waves (E = 0), the N = 1 state has the lowest value of the Lyapunov 
function and hence will be the selected pattern. At low frequency, the system effectively 
crosses over to the large damping region regardless of its (finite) viscosity (this range 
was not accessible to the low damping calculation of of Zhang & Vinals (1997)). On the 
other hand, for pure capillary waves (E = 1) the preferred pattern is N — 2 at low 
damping and TV = 1 at high damping. Interesting behavior is observed in the vicinity 
of E = 1/3 (mixed gravity capillary waves) where the triad resonance angle approaches 
zero. Hexagonal and higher symmetry quasipatterns are selected with decreasing 7. The 
low damping results in this region are in qualitative agreement with the earlier work of 
Zhang fc Vinals (1997) , although the latter could not account for the transition to N = 1 
as 7 is increased. 



We finally compare our predictions (based on Eq. (4.1)) and two recent sets of experi- 



ments that addressed pattern selection in the large aspect ratio limit by Kudrolli 
and by 



Gollub (1996) 



Binks fc van de Water (1997) . The only input parameters in our calculations are 



the fluid properties (density, surface tension and viscosity), and the frequency of the 
driving acceleration. All these parameters are known fairly precisely in the experimental 
work, thus allowing a quantitative comparison between theory and experiments. 



Pattern selection in the low viscosity range has been recently studied by Binks & van de Water (1997) 
They have developed a cell of exceptionally large aspect ratio, and of depth that is much 
larger than the wavelength. The fluid used was a low viscosity, low surface tension silicon 
oil with v = 0.03397 cm 2 /s, p — 0.8924 g/cm 3 and a — 18.3 dyne/cm. Given the range of 
frequencies studied, the viscous damping parameter probed was within 7 ~ 0.01 — 0.03. 
They have reported transitions from standing wave patterns of square symmetry at high 
frequency (l 41 Hz), to hexagonal, eight-fold and ten- fold quasi-periodic patterns upon 
lowering the driving frequency. Stable hexagonal patterns appear at approximately 
36 Hz, although a transition region of mixed square/hexagonal symmetry is observed 
between approximately 36 Hz and 41 Hz. Given the parameters of this experiment, our 
theory predicts a transition at 35.4 Hz , compared to the value of 32.8 Hz given by the 
earlier work of Zhang fc Vinals (1997) . An additional transition region exhibiting pat- 



terns of mixed hexagonal and eight-fold symmetry was also observed between 30-31 Hz, 
which compares favorab ly with our prediction for the transition to eight-fold symmet- 
ric patterns at 28.7 Hz ( Zhang fc Vinals (1997) had predicted the transition to occur at 
27.9 Hz). [j] A possible explanation for the lar ger discrepancy between the experiments 
and the calculations of Zhang fc Vinals (1997) involves the fact that t hey o nly used the 
term linear in 7 in the calculation of the threshold for instability (Eq. ( 2.14 )). Omitting 
the first correction (of order 7 3 ' 2 ) yields a similar percentage error in the threshold value 
(for 7 in the range 0.01 - 0.03 as is appropriate for this experiment). 

Another set of recent experiments involving fluids of different viscosity has been carried 
out by Kudrolli fc Gollub (1996)| . Although the depth of the fluid layer (0.3 cm) is smaller 
than the wavelength of the waves (1-3 cm), the comparison is still illuminating. Figure 
H shows the symmetry of the predicted patterns as a function of the viscosity of the 
fluid and of the driving frequency (with p — 0.95g/cm 3 and a = 20.6dyn/cm), and 
the experimentally observed patterns. They find a stripe pattern at high viscosity, a 
hexagonal pattern at low viscosity and frequency, and a square pattern at low viscosity 
and high frequency. Two significant discrepancies concern the experimental observation 
of a hexagonal pattern at v = 1cm 2 /s and low frequency, and also at v = 0.04cm 2 /s 



f In this experiment, patterns with N = 5 are observed around 27Hz. We agree with the 
authors of the experiment that this discrepancy may be due to finite size effects. In fact, J-{5) 
is very close to .F(4) at about 26Hz (the difference is less than 0.2%), although J- (5) > ^(4). 
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and / = 27Hz. It is possible that the shallowness of the fluid layer can account for 
these differences, especially in view of the fact that, as noted above, the experiments by 



Binks & van de Water (1997) did probe this latter region in a deep fluid layer, and their 
results do agree with our predictions. 

Finally, the fact that portions of the boundaries separating regions of different sym- 
metry appear almost as straight lines in Fig. pi is due to the log-log scale used. In 
addition, the transition line delimiting the region of stripe patterns is almost indepen- 
dent of frequency only because of the limited frequency range probed in the experiments 
and displayed in the figure. On the other hand, the line separating regions of square 
and hexagonal patterns is almost independent of viscosity because it depends mainly on 
whether the waves are capillarity or gravity dominated, fact that is largely dependent on 
the driving frequency and not on viscosity. 

In summary, we have presented a nonlinear theory for Faraday waves in viscous fluids 
with no assumptions or approximations other than those inherent to the multi-scale 
expansion. A set of standing wave amplitude equations has been obtained that is of 
gradient form. Minimization of the associated Lyapunov function leads to determination 
of the preferred pattern near threshold. The predicted patterns are in excellent agreement 
with recent experiments in large aspect ratio systems involving a range of fluid viscosities 
and driving frequencies. According to Fig. n], the transition from square to stripe patterns 
remains in the capillary wave limit of S = 1 (high frequency limit in the experiments). 
However, the figure for £ = indicates that stripe patterns are always preferred in the 
pure-gravity- wave limit (low frequency limit in the experiments). Furthermore, all the 
high symmetry patterns (with N > 3) are observed in the vicinity of £ = 1/3, point at 
which the triad resonant angle approaches zero, and for low damping where the resonance 
is more pronounced. 

This research has been supported by the U.S. Department of Energy, contract No. 
DE-FG05-95ER14566, and also in part by the Supercomputer Computations Research 
Institute, which is partially funded by the U.S. Department of Energy, contract No. 
DE-FC05-85ER25000. 
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Appendix A. Inhomogeneous terras of the first and second order 
equations 

We list in this appendix the functions Gij , the inhomogeneous terms in the boundary 
conditions at first and second order. 

Gn = d z woCa - u d x Co - v d y (o, 

G12 = d x [-d zz u (o - d xz w a ( + 2(d x u - d z w )d x ( a + (d y u + d x v )dy( ] 

+ d y [-d zz v Co - d yz w (o + 2(d y v - d z w a )d v Co + (d x v + d y u Q )d x ( } , 
G13 = -pVij • [(uo • V)u ] + \7 2 H (-2nd zz w (o + d z p (o) 
G 2 \ = -d T (o - uod x Ci - uid x Co - d z u (od x Co - v a d y C,i - vid y ( a - d z v C, d y C, 

+ d z w (i + d z wx(o + \d zz w C,l, 



G22 — d x 



- d zz uiCo - d zz u Ci - hd zzz u (Q - d xz wi( Q - d xz w (i - ±d xzz w (;% 



- 2(d z w 1 - d x ui)d x ( - 2(d z w - d x u )d x (i - 2d z (d z w - d x u )C d x ( 
+ (d y ui + d x vi)d y Co + (d y u + d x v )d y Ci + d z (d y u a + d x v a )( d y C 

+ d y - d^uiCo - <9 zz v Ci ~ ^d zzz v Co - dyzWiCo ~ d yz w Ci - \d yzz w C,l 

- 2(d z wi ~ d y vi)dyCo - 2(d z w - d y v )d y (i - 2d z (d z w - d y v )Cod y Co 
+ (d y ui + d x vi)d x ( + {dyu + d x v )d x C,i + d z (d v u + d x v )( d x (o 



G 23 = -pV ff • [(uo • V)ui + (ui • Vu )] + pd T d z w + V 2 H [ - \pf a {e 11 
+ d z piC, + d z p d + ^d z p (o - 2 vd z W!(o - 2rjd z w (i - vd z w ( 2 
+ 2i](d z ui + d x wi)d x ( + 2r)(d z w - d x u )(d x ( ) 2 
+ 2i](d z v 1 + dytu^dyCo + 2r](d z w - d y v )(d y ( ) 2 

- 2i](d y u + d x v )d x ( dyCo ~ ^d xx Co{d x ( ) 2 - %vd yy (o{dy( ) 2 

- \vd xx C,o{d y (o) 2 - ±ad yy ( (d x (o) 2 - 2ad x (odyCod xv (o\. 



e-™) Co 
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Appendix B. Matrix of coefficients at first order 

Left hand side of the system of linear equations for the first order solution (for sim- 
plicity, we only show the case 7 ^C 1 and the coefficients of first time harmonic e luJt / 2 ), 



( -1 








-1 








-1% 





^ 




fa 1 - \ 





-1 








-1 
















a 








-1 








-1 








2i 




a 1+ 


- 7 fc 2 
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l ~ 





-p 








-k 2 
















0° 
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-2i 













1+ 
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-2ik 








7fc 2 Q* 
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27A; 2 
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^mn 
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Gk 2 + S/c 4 


2 7 /c 2 




6° 

mn 


I 





2ik 








jk 2 q 
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Figure 1 . Schematic setup of a Faraday wave configuration. 
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Figure 2. Dimensionless threshold for linear instability A on sct as a function of the dimensionless 
damping parameter 7. The lower solid line is the exact result; the upper solid line is the lowest 
order approximation in the damping parameter. Also shown are the first order correction in the 
viscous damping parameter (dashed line), and the first correction for the instability threshold 
for a damped Mathieu equation. 
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Figure 3. Cubic term coefficient of the standing wave amplitude equation as a function of 
angle between wavevectors 6, in the limit of gravity waves, E = 0, and different viscous damping 
coefficients. 
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Figure 4. Cubic term coefficient of the standing wave amplitude equation as a function of 
angle between wavevectors 6, in the mixed capillary-gravity regime (E = 1/3). Note that the 
curve becomes extremely flat near cos# = for low 7. 
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Figure 5. Schematic representation of a triad resonant interaction: two linearly unstable modes 
k m and k n interact to produce a linearly stable mode. This mode interacts with — k n leading 
to resonance with k m . 
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Figure 6. Cubic term coefficient of the standing wave amplitude equation as a function of 
angle between wavevectors 6, in the limit of capillary waves, E = 1. The large peaks at small 
values of 7 are due to triad resonant interactions. 
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Figure 7. Numerical values of the Lyapunov function for regular patterns comprising N stand- 
ing waves as a function of the viscous damping parameter 7. Bottom right: gravity wave limit; 
bottom left: capillary wave limit; top left, the mixed case of £ = 1/3; top right is the same as 
top left but showing the region of small damping in more detail. In the two bottom plots, the 
curves not labeled are ordered in increasing order of N. 
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Figure 8. Selected patterns a s a function of fluid viscos ity and driving frequency. The symbols 
are the experimental results of Kudrolli & Gollub (1996). x, stripe patterns; □, square patterns; 
and, A, hexagonal patterns. Alternating x and □ indicate regions in which stationary mixed 
stripe and square patterns were observed. 



